# 文件创建时间: 2025-09-24
# 文件功能描述: 稳态两流体求解器，实现Picard迭代求解两相流动量和压力修正方程

import numpy as np
from typing import Dict, List, Tuple
import os
import logging
from physics import PhaseType
from flow_regimes import FlowRegimeIdentifier
from multigrid import MultiGridSolver
from scipy.sparse import lil_matrix, bmat, csr_matrix, diags
from scipy.linalg import solve_banded
from scipy.sparse.linalg import bicgstab, spsolve, gmres
import scipy.sparse
from collections import deque

class PhysicsBasedPreconditioner:
    """基于物理的预处理器"""
    def __init__(self):
        self.matrix = None
        
    def build_matrix(self, system_matrix, subchannel):
        """构建预处理器矩阵
        
        Args:
            system_matrix: 系统矩阵
            subchannel: 子通道对象
            
        Returns:
            np.ndarray: 预处理器矩阵
        """
        n = system_matrix.shape[0]
        self.matrix = np.zeros_like(system_matrix)
        
        # 提取主对角线及其相邻元素
        for i in range(n):
            # 主对角线
            self.matrix[i,i] = system_matrix[i,i]
            
            # 下对角线
            if i > 0:
                self.matrix[i,i-1] = system_matrix[i,i-1]
            
            # 上对角线
            if i < n-1:
                self.matrix[i,i+1] = system_matrix[i,i+1]
        
        # 确保矩阵非奇异
        for i in range(n):
            self.matrix[i,i] = max(self.matrix[i,i], 
                                 1.1 * np.sum(np.abs(self.matrix[i,:])))
        
        return self.matrix

class SteadyStateSolver:
    """稳态两流体求解器"""
    
    def __init__(self, grid, physics):
        """初始化求解器
        
        Args:
            grid: 计算网格
            physics: 物理模型
        """
        self.grid = grid
        self.physics = physics
        self.boundary_conditions = {}
        
        # Picard迭代参数
        self.max_iterations = 1000
        self.convergence_tolerance = 1e-4
        
        # 阶段3优化：添加压力修正频率控制（每2-3次动量迭代执行一次压力修正）
        self.pressure_correction_frequency = 1 # 改为1，每次都修正
        
        # 伪瞬态时间步长 (CFL数)
        self.pseudo_cfl = 0.5
        # NOTE: a_P 缓存机制，为 PKM简化算法提供速度修正系数
        self.momentum_ap = {}
        
        # 初始化松弛因子
        self._setup_relaxation_factors()
        
        # 物理限制
        self.min_pressure = 1e5      # 最小压力 (1 bar)
        self.max_pressure = 20e6     # 最大压力 (200 bar)
        self.min_velocity = -20.0    # 最小速度
        self.max_velocity = 20.0     # 最大速度
        self.min_void = 1e-6         # 最小空泡份额
        self.max_void = 1 - 1e-6     # 最大空泡份额
        
        # GMRES求解器参数
        self.gmres_tolerance = 1e-6   # 提高GMRES精度
        self.max_gmres_iterations = 500  # 增加最大迭代次数
        self.gmres_restart = 100     # 增加重启参数
        
        # 重力加速度
        self.g = 9.81  # m/s²
        
        # 初始化预处理器
        self.preconditioner = PhysicsBasedPreconditioner()
        
        # 使用物理模型中的流型识别器
        self.flow_regime_identifier = self.physics.flow_regime_identifier
        
        # 为SIMPLE算法引入中间速度场
        self.intermediate_velocity = None

    def _setup_relaxation_factors(self):
        """设置自适应松弛因子策略
        
        【修改2025-11-05】支持分离的气液相速度松弛因子
        - velocity_g: 气相速度松弛因子（默认0.2，更保守以提高稳定性）
        - velocity_l: 液相速度松弛因子（默认0.4，较大值以保证收敛速度）
        """
        # 【优化】调整基础松弛因子，平衡收敛速度和稳定性
        self.base_relaxation_factors = {
            'pressure': 0.1,           # 适度提高压力松弛因子（原0.05太保守）
            'velocity_g': 0.2,         # 【新增】气相速度松弛因子（参考COBRA的DAMPFF机制）
            'velocity_l': 0.4,         # 【新增】液相速度松弛因子
            'void_fraction': 0.15,     # 提高空泡份额松弛因子（原0.1）
            'velocity_correction': 0.2 # 提高速度修正松弛因子（原0.1）
        }
        
        # 当前使用的松弛因子
        self.relaxation_factors = self.base_relaxation_factors.copy()
        
        # TODO: 任务 2.3 - 使用deque限制残差历史记录长度，防止内存泄漏
        # 残差历史记录
        self.residual_history = {
            'pressure': deque(maxlen=20),
            'velocity_g': deque(maxlen=20),
            'velocity_l': deque(maxlen=20),
            'mass_g': deque(maxlen=20),
            'mass_l': deque(maxlen=20)
        }
        
        # 残差振荡检测阈值
        self.oscillation_threshold = 1.2  # 如果残差比前一步高出20%，认为有振荡
        
        # 松弛因子调整参数
        self.relax_decrease_factor = 0.8  # 松弛因子减小倍数
        self.relax_increase_factor = 1.05  # 松弛因子增加倍数
        self.min_relaxation = 0.01  # 最小松弛因子
        self.max_relaxation = 0.5   # 最大松弛因子

    def _adjust_relaxation_factors(self, residuals):
        """根据残差历史动态调整松弛因子
        
        【修复 2025-11-05】临时禁用自适应调整，避免假收敛场景下的误判
        
        Args:
            residuals: 包含各类残差的字典，键为'pressure', 'velocity_g', 'velocity_l', 'mass_g', 'mass_l'
        """
        # 临时禁用自适应调整，直到假收敛问题解决
        return
        # 更新残差历史
        self.residual_history['pressure'].append(residuals['pressure'])
        self.residual_history['velocity_g'].append(residuals['velocity_g'])
        self.residual_history['velocity_l'].append(residuals['velocity_l'])
        self.residual_history['mass_g'].append(residuals['mass_g'])
        self.residual_history['mass_l'].append(residuals['mass_l'])
        
        # 至少需要2步历史才能调整
        if len(self.residual_history['pressure']) < 2:
            return
        
        # 检查各残差是否振荡
        # 【修复】只有在残差足够大时才检测振荡，避免在接近收敛时误判
        oscillations = {}
        for key in self.residual_history:
            current = self.residual_history[key][-1]
            previous = self.residual_history[key][-2]
            
            # 只有当残差大于收敛标准时才检测振荡
            # 避免在残差已经很小时，微小波动被误判为振荡
            threshold_map = {
                'pressure': 1e-2,
                'velocity_g': 1e-3,
                'velocity_l': 1e-3,
                'mass_g': 1e-3,
                'mass_l': 1e-2
            }
            
            if previous < threshold_map.get(key, 1e-3):
                # 残差已经很小，不检测振荡
                oscillations[key] = False
            else:
                ratio = current / previous
                oscillations[key] = ratio > self.oscillation_threshold
        
        # 根据振荡情况调整松弛因子
        if oscillations['pressure']:
            self.relaxation_factors['pressure'] = max(
                self.relaxation_factors['pressure'] * self.relax_decrease_factor, 
                self.min_relaxation
            )
            logging.info(f"检测到压力残差振荡，降低压力松弛因子至{self.relaxation_factors['pressure']:.4f}")
        
        # 【修改2025-11-05】分别调整气相和液相速度松弛因子
        if oscillations['velocity_g']:
            self.relaxation_factors['velocity_g'] = max(
                self.relaxation_factors['velocity_g'] * self.relax_decrease_factor, 
                self.min_relaxation
            )
            logging.info(f"检测到气相速度残差振荡，降低气相速度松弛因子至{self.relaxation_factors['velocity_g']:.4f}")
        
        if oscillations['velocity_l']:
            self.relaxation_factors['velocity_l'] = max(
                self.relaxation_factors['velocity_l'] * self.relax_decrease_factor, 
                self.min_relaxation
            )
            logging.info(f"检测到液相速度残差振荡，降低液相速度松弛因子至{self.relaxation_factors['velocity_l']:.4f}")
        
        if oscillations['mass_g'] or oscillations['mass_l']:
            self.relaxation_factors['void_fraction'] = max(
                self.relaxation_factors['void_fraction'] * self.relax_decrease_factor, 
                self.min_relaxation
            )
            logging.info(f"检测到质量残差振荡，降低空泡份额松弛因子至{self.relaxation_factors['void_fraction']:.4f}")
        
        # 如果残差持续减小，可以尝试增加松弛因子（加速收敛）
        for key in self.residual_history:
            if len(self.residual_history[key]) >= 4:  # 需要更多历史判断趋势
                last_3_decreasing = True
                for i in range(1, 3):
                    if self.residual_history[key][-i-1] <= self.residual_history[key][-i]:
                        last_3_decreasing = False
                        break
                
                if last_3_decreasing:
                    if key == 'pressure':
                        self.relaxation_factors['pressure'] = min(
                            self.relaxation_factors['pressure'] * self.relax_increase_factor, 
                            self.max_relaxation
                        )
                    elif key == 'velocity_g':
                        self.relaxation_factors['velocity_g'] = min(
                            self.relaxation_factors['velocity_g'] * self.relax_increase_factor, 
                            self.max_relaxation
                        )
                    elif key == 'velocity_l':
                        self.relaxation_factors['velocity_l'] = min(
                            self.relaxation_factors['velocity_l'] * self.relax_increase_factor, 
                            self.max_relaxation
                        )
                    elif key in ['mass_g', 'mass_l']:
                        self.relaxation_factors['void_fraction'] = min(
                            self.relaxation_factors['void_fraction'] * self.relax_increase_factor, 
                            self.max_relaxation
                        )
    
    def _apply_numerical_stabilization(self, subchannels):
        """应用增强的数值稳定化技术"""
        logging.info("应用增强的数值稳定化技术")
        
        max_alpha_change = 0.05
        max_pressure_rel_change = 0.02
        
        max_velocity_change = {
            'vapor': {'axial': 0.3, 'lateral': 0.2},
            'liquid': {'axial': 0.2, 'lateral': 0.1}
        }
        
        for subchannel in subchannels:
            delta_alpha = subchannel.alpha - subchannel.previous_alpha
            limited_delta_alpha = np.clip(delta_alpha, -max_alpha_change, max_alpha_change)
            subchannel.alpha = subchannel.previous_alpha + limited_delta_alpha
            
            delta_p_rel = (subchannel.pressure - subchannel.previous_pressure) / np.maximum(subchannel.previous_pressure, 1.0)
            limited_delta_p_rel = np.clip(delta_p_rel, -max_pressure_rel_change, max_pressure_rel_change)
            subchannel.pressure = subchannel.previous_pressure * (1.0 + limited_delta_p_rel)
            
            for phase in ['vapor', 'liquid']:
                for direction in ['axial', 'lateral']:
                    v_curr = subchannel.velocity[phase][direction]
                    v_prev = subchannel.previous_velocity[phase][direction]
                    delta_v = v_curr - v_prev
                    max_change = max_velocity_change[phase][direction]
                    
                    # 【修复】只在第一个通道打印日志，避免重复
                    if phase == 'liquid' and direction == 'axial' and hasattr(self, 'residual_history') and len(self.residual_history['mass_l']) > 2:
                        if self.residual_history['mass_l'][-1] > self.residual_history['mass_l'][-2]:
                            max_change *= 0.5
                            if subchannel.id == 0:
                                logging.info("检测到液相质量残差增加，增强液相速度稳定化")
                    
                    limited_delta_v = np.clip(delta_v, -max_change, max_change)
                    subchannel.velocity[phase][direction] = v_prev + limited_delta_v
            
            # TODO: 任务 2.3 - 改进数值稳定化措施，对alpha进行更严格的物理限制
            subchannel.alpha = np.clip(subchannel.alpha, self.min_void, 1.0 - self.min_void)
            subchannel.pressure = np.clip(subchannel.pressure, self.min_pressure, self.max_pressure)
            
            for phase in ['vapor', 'liquid']:
                for direction in ['axial', 'lateral']:
                    subchannel.velocity[phase][direction] = np.clip(
                        subchannel.velocity[phase][direction], 
                        self.min_velocity, self.max_velocity)
        
        self._log_field_statistics(subchannels)
    
    def _log_field_statistics(self, subchannels):
        """记录场量统计信息，用于监控稳定性"""
        min_alpha, max_alpha = float('inf'), float('-inf')
        min_pressure, max_pressure = float('inf'), float('-inf')
        
        velocity_stats = {
            'vapor': {'axial': {'min': float('inf'), 'max': float('-inf')}, 'lateral': {'min': float('inf'), 'max': float('-inf')}},
            'liquid': {'axial': {'min': float('inf'), 'max': float('-inf')}, 'lateral': {'min': float('inf'), 'max': float('-inf')}}
        }
        
        for subchannel in subchannels:
            min_alpha = min(min_alpha, np.min(subchannel.alpha))
            max_alpha = max(max_alpha, np.max(subchannel.alpha))
            min_pressure = min(min_pressure, np.min(subchannel.pressure))
            max_pressure = max(max_pressure, np.max(subchannel.pressure))
            
            for phase in ['vapor', 'liquid']:
                for direction in ['axial', 'lateral']:
                    v = subchannel.velocity[phase][direction]
                    velocity_stats[phase][direction]['min'] = min(velocity_stats[phase][direction]['min'], np.min(v))
                    velocity_stats[phase][direction]['max'] = max(velocity_stats[phase][direction]['max'], np.max(v))
        
        logging.info(f"场量统计 - 空泡份额: [{min_alpha:.4f}, {max_alpha:.4f}]")
        logging.info(f"场量统计 - 压力: [{min_pressure:.2e}, {max_pressure:.2e}] Pa")
        logging.info(f"场量统计 - 气相轴向速度: [{velocity_stats['vapor']['axial']['min']:.2f}, {velocity_stats['vapor']['axial']['max']:.2f}] m/s")
        logging.info(f"场量统计 - 液相轴向速度: [{velocity_stats['liquid']['axial']['min']:.2f}, {velocity_stats['liquid']['axial']['max']:.2f}] m/s")
            
    def solve_steady_state(self, subchannels):
        """求解稳态问题"""
        logging.info("\n开始稳态求解")
        self._save_current_fields(subchannels)
        
        max_residual = float('inf')
        for iter_num in range(self.max_iterations):
            logging.info(f"\n第 {iter_num + 1} 次迭代:")
            
            res_ax_g, res_ax_l, res_lat_g, res_lat_l = self._solve_steady_momentum(subchannels)
            logging.info(f"  气相轴向动量残差: {res_ax_g:.3e}")
            logging.info(f"  液相轴向动量残差: {res_ax_l:.3e}")
            logging.info(f"  气相横向动量残差: {res_lat_g:.3e}")
            logging.info(f"  液相横向动量残差: {res_lat_l:.3e}")
            
            # 3. [PKM简化] 耦合求解压力和空泡份额
            # 采用降维+Krylov方法处理气相和液相的连续性方程，并更新压力、空泡份额和速度场
            pressure_residual, void_fraction_residual, success = self._solve_pressure_void_coupled(subchannels)
            if not success:
                logging.error("耦合求解器失败，中止计算。")
                return False

            # 4. 使用真实的质量守恒残差
            # FIXME: 任务 1.2 - 使用真实的质量残差
            res_mass_g = np.max(self.mass_residual_g) if hasattr(self, 'mass_residual_g') else 0.0
            res_mass_l = np.max(self.mass_residual_l) if hasattr(self, 'mass_residual_l') else 0.0
            logging.info(f"  气相质量守恒残差: {res_mass_g:.3e}")
            logging.info(f"  液相质量守恒残差: {res_mass_l:.3e}")
            
            self._apply_numerical_stabilization(subchannels)
            
            # 【修改2025-11-05】实现分层收敛判据，参考COBRA的实践
            # 不同物理量使用不同的收敛标准
            velocity_g_converged = res_ax_g < 5e-3  # 气相速度放宽到5e-3（考虑固有振荡）
            velocity_l_converged = res_ax_l < 1e-3  # 液相速度保持严格标准
            lateral_g_converged = res_lat_g < 1e-2  # 横向动量放宽标准（通常难收敛）
            lateral_l_converged = res_lat_l < 5e-3
            pressure_converged = pressure_residual < 1e-2
            void_fraction_converged = void_fraction_residual < 1e-4
            mass_g_converged = res_mass_g < 1e-3
            mass_l_converged = res_mass_l < 1e-3
            
            # 综合判断
            all_converged = (velocity_g_converged and velocity_l_converged and 
                           lateral_g_converged and lateral_l_converged and
                           pressure_converged and void_fraction_converged and 
                           mass_g_converged and mass_l_converged)
            
            # 保留最大残差用于日志输出
            max_residual = max(res_ax_g, res_ax_l, res_lat_g, res_lat_l, 
                              pressure_residual, void_fraction_residual, res_mass_g, res_mass_l)
            
            if all_converged:
                logging.info(f"\n求解在第 {iter_num + 1} 次迭代收敛")
                logging.info(f"  气相轴向速度: {res_ax_g:.3e} < 5e-3 ✓")
                logging.info(f"  液相轴向速度: {res_ax_l:.3e} < 1e-3 ✓")
                logging.info(f"  气相横向动量: {res_lat_g:.3e} < 1e-2 ✓")
                logging.info(f"  液相横向动量: {res_lat_l:.3e} < 5e-3 ✓")
                logging.info(f"  压力: {pressure_residual:.3e} < 1e-2 ✓")
                logging.info(f"  空泡份额: {void_fraction_residual:.3e} < 1e-4 ✓")
                logging.info(f"  气相质量: {res_mass_g:.3e} < 1e-3 ✓")
                logging.info(f"  液相质量: {res_mass_l:.3e} < 1e-3 ✓")
                logging.info(f"最大残差: {max_residual:.3e}")
                return True
                
            self._save_current_fields(subchannels)
            
            residuals = {
                'pressure': pressure_residual,
                'velocity_g': res_ax_g,
                'velocity_l': res_ax_l,
                'mass_g': res_mass_g,
                'mass_l': res_mass_l
            }
            self._adjust_relaxation_factors(residuals)
            
        logging.warning(f"\n求解未在 {self.max_iterations} 次迭代内收敛")
        logging.warning(f"最终残差: {max_residual:.3e}")
        return False
        
    def _save_current_fields(self, subchannels):
        """保存当前场值用于计算残差"""
        for subchannel in subchannels:
            subchannel.previous_pressure = np.copy(subchannel.pressure)
            subchannel.previous_alpha = np.copy(subchannel.alpha)
            subchannel.previous_mass_flux = {
                'vapor': np.zeros_like(subchannel.pressure),
                'liquid': np.zeros_like(subchannel.pressure)
            }
            subchannel.previous_velocity = {
                'vapor': {'axial': np.copy(subchannel.velocity['vapor']['axial']), 'lateral': np.copy(subchannel.velocity['vapor']['lateral'])},
                'liquid': {'axial': np.copy(subchannel.velocity['liquid']['axial']), 'lateral': np.copy(subchannel.velocity['liquid']['lateral'])}
            }
            
    def _solve_steady_momentum(self, subchannels):
        """求解稳态动量方程，并将结果存入中间速度场"""
        # 初始化中间速度场
        self.intermediate_velocity = {
            ch.id: {
                'vapor': {'axial': np.copy(ch.velocity['vapor']['axial']), 'lateral': np.copy(ch.velocity['vapor']['lateral'])},
                'liquid': {'axial': np.copy(ch.velocity['liquid']['axial']), 'lateral': np.copy(ch.velocity['liquid']['lateral'])}
            } for ch in subchannels
        }

        max_residual_axial_g, max_residual_axial_l = 0.0, 0.0
        max_residual_lateral_g, max_residual_lateral_l = 0.0, 0.0
        n_axial = self.grid.n_axial
        
        # 添加详细的调试信息
        logging.info("  开始求解动量方程")
        
        # 检查输入数据的有效性
        nan_count = 0
        inf_count = 0
        for subchannel in subchannels:
            if np.any(np.isnan(subchannel.pressure)) or np.any(np.isinf(subchannel.pressure)):
                nan_count += 1
            if np.any(np.isnan(subchannel.alpha)) or np.any(np.isinf(subchannel.alpha)):
                inf_count += 1
        
        if nan_count > 0 or inf_count > 0:
            logging.warning(f"  输入数据异常: {nan_count}个通道有NaN/Inf压力, {inf_count}个通道有NaN/Inf空泡份额")

        for phase in ['vapor', 'liquid']:
            for subchannel in subchannels:
                diag, lower, upper, rhs = np.zeros(n_axial), np.zeros(n_axial-1), np.zeros(n_axial-1), np.zeros(n_axial)
                
                # 详细记录第一个通道的计算过程
                debug_channel = (subchannel.id == 0)
                if debug_channel:
                    logging.info(f"    调试通道{subchannel.id} {phase}相动量方程:")
                
                v_old = subchannel.velocity[phase]['axial'].copy()

                for i in range(n_axial):
                    # 首先获取当前节点的物性和空泡份额
                    props = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.VAPOR if phase == 'vapor' else PhaseType.LIQUID)
                    rho = props['density']
                    mu = props['viscosity']
                    alpha = subchannel.alpha[i]

                    # --- 时间项 (参考COBRA: DXDT=DX/DT) ---
                    # 【修复 2025-10-30v3】量纲统一：v_new/dt = v_old/dt + dt*acceleration
                    # DIAG和RHS的量纲都是速度 m/s
                    alpha_phase = alpha if phase == 'vapor' else (1 - alpha)
                    # 【修复 2025-11-05】进一步降低气相伪时间步长以提高稳定性
                    pseudo_dt = 1e-4 if phase == 'vapor' else 1e-3  # 气相1e-4（从5e-4进一步降低），液相1e-3
                    # 时间项系数（量纲：1/s）
                    time_coeff = 1.0 / pseudo_dt
                    diag[i] += time_coeff
                    rhs[i] += time_coeff * v_old[i]

                    # Convection term
                    if i > 0:
                        v_up = subchannel.velocity[phase]['axial'][i-1]
                        if v_up > 0:
                            lower[i-1] = -v_up
                            diag[i] += v_up
                        else:
                            diag[i-1] += v_up
                            upper[i-1] = -v_up
                    
                    # 检查物性计算结果
                    if debug_channel and i < 3:
                        logging.info(f"      节点{i}: P={subchannel.pressure[i]:.2e}Pa, alpha={alpha:.4f}, rho={rho:.2f}, mu={mu:.2e}")
                    
                    if not (np.isfinite(rho) and np.isfinite(mu) and np.isfinite(alpha)):
                        logging.error(f"    通道{subchannel.id}节点{i}物性异常: rho={rho}, mu={mu}, alpha={alpha}")
                        continue
                    
                    # Pressure term - 【修复】乘以dt使量纲为m/s
                    if i < n_axial - 1 and rho > 1e-6:
                        pressure_gradient = (subchannel.pressure[i+1] - subchannel.pressure[i]) / self.grid.dz
                        # 压力项量纲：(m/s²) * s = m/s ✓
                        rhs[i] -= pseudo_dt * (1/rho) * pressure_gradient
                        
                        if debug_channel and i < 3:
                            logging.info(f"      节点{i}: 压力梯度项={pressure_gradient:.2e}, 贡献={-pseudo_dt*pressure_gradient/rho:.2e}")
                    
                    # Gravity term - 【修复】乘以dt使量纲为m/s
                    # 量纲：(m/s²) * s = m/s ✓
                    rhs[i] -= pseudo_dt * self.g * self.grid.cos_theta
                    
                    # 【修复 2025-10-30】参考COBRA：摩擦项隐式处理（放在DIAG）
                    # COBRA: DPKL = f*DX*ρ*v/D_h*|v|/A，放在DIAG中
                    # 这样摩擦力会自动平衡加速，提供数值稳定性
                    v = subchannel.velocity[phase]['axial'][i]
                    
                    # 计算摩擦系数（隐式形式）
                    if abs(v) > 1e-6 and rho > 1e-6:
                        # 摩擦系数 = f * P_w / (2 * A * D_h) * |v|
                        # 其中f可以从雷诺数估算
                        Re = rho * abs(v) * subchannel.hydraulic_diameter / (mu + 1e-10)
                        if Re > 2300:  # 湍流
                            f = 0.316 / (Re**0.25)  # Blasius公式
                        else:  # 层流
                            f = 64.0 / (Re + 1e-6)
                        
                        # 隐式摩擦项系数（放在DIAG）
                        # 摩擦力造成的减速度：f * P_w * |v| / (2*A*D_h*ρ)
                        # 乘以dt使量纲为 1/s
                        friction_coeff = pseudo_dt * 0.5 * f * subchannel.wetted_perimeter * abs(v) / (subchannel.area * subchannel.hydraulic_diameter)
                        diag[i] += friction_coeff
                        
                        if debug_channel and i < 3:
                            logging.info(f"      节点{i}: 摩擦系数={friction_coeff:.2e}, Re={Re:.0f}, f={f:.4f}")

                    # 【修复 2025-10-30】界面摩擦项也隐式处理
                    # COBRA方法：相间摩擦力与相对速度成正比，隐式处理提供稳定性
                    v_other = subchannel.velocity['liquid' if phase == 'vapor' else 'vapor']['axial'][i]
                    
                    # 简化界面摩擦模型：F_i ∝ (v_g - v_l)
                    # 隐式形式：将F_i/v项放在DIAG，F_i中的v_other项放在RHS
                    if alpha > 1e-6 and (1-alpha) > 1e-6:
                        # 界面摩擦系数 - 【修复 2025-11-04】使用物理模型中的值
                        # 从物理模型中获取界面摩擦系数，而非硬编码
                        C_i = self.physics.C_i  # 界面摩擦系数
                        interfacial_coeff = C_i * subchannel.wetted_perimeter * self.grid.dz
                        
                        # 【修复】界面摩擦也乘以dt统一量纲
                        interfacial_coeff_dt = pseudo_dt * interfacial_coeff
                        if phase == 'vapor':
                            # 气相：F_i = C_i * (v_g - v_l)
                            # 隐式：diag += dt*C_i, rhs -= dt*C_i * v_l
                            diag[i] += interfacial_coeff_dt / (alpha * rho + 1e-10)
                            rhs[i] -= interfacial_coeff_dt * v_other / (alpha * rho + 1e-10)
                        else:
                            # 液相：F_i = -C_i * (v_g - v_l)
                            diag[i] += interfacial_coeff_dt / ((1-alpha) * rho + 1e-10)
                            rhs[i] += interfacial_coeff_dt * v_other / ((1-alpha) * rho + 1e-10)
                        
                        if debug_channel and i < 3:
                            logging.info(f"      节点{i}: 界面摩擦系数={interfacial_coeff:.2e}")

                # 【修复 2025-10-30】移除人工数值扩散
                # COBRA不使用人工扩散，而是依靠物理摩擦项提供稳定性
                # 人工扩散会破坏物理精度，导致非物理的速度分布
                # nu = 0.2  # 已移除
                # diag += 2 * nu / (self.grid.dz ** 2)
                # if n_axial > 1:
                #     upper -= nu / (self.grid.dz ** 2)
                #     lower -= nu / (self.grid.dz ** 2)
                
                # 检查系数矩阵的合理性
                if debug_channel:
                    logging.info(f"    系数矩阵诊断: diag范围=[{np.min(diag):.2e}, {np.max(diag):.2e}]")
                    logging.info(f"    右端项诊断: rhs范围=[{np.min(rhs):.2e}, {np.max(rhs):.2e}]")
                    # 【新增】检查对角占优性
                    if n_axial > 2:
                        diag_dominance = np.min(np.abs(diag[1:-1]) / (np.abs(upper[:-1]) + np.abs(lower[1:]) + 1e-10))
                        logging.info(f"    对角占优度: {diag_dominance:.2e} (>1为对角占优)")
                
                # 检查矩阵条件数
                if np.any(np.abs(diag) < 1e-12):
                    logging.warning(f"    通道{subchannel.id} {phase}相: 系数矩阵对角元过小")
                    diag = np.maximum(np.abs(diag), 1e-10) * np.sign(diag)
                
                ab = np.zeros((3, n_axial)); ab[0, 1:] = upper; ab[1, :] = diag; ab[2, :-1] = lower
                
                # NOTE: 缓存 a_P 系数，供 PKM简化耦合求解器使用
                if subchannel.id not in self.momentum_ap:
                    self.momentum_ap[subchannel.id] = {'vapor': {}, 'liquid': {}}
                self.momentum_ap[subchannel.id][phase]['axial'] = diag.copy()

                # 【修改2025-11-05】对气相和液相使用不同的松弛因子
                omega = self.relaxation_factors['velocity_g'] if phase == 'vapor' else self.relaxation_factors['velocity_l']
                
                try:
                    v_new = solve_banded((1, 1), ab, rhs)
                except Exception as e:
                    logging.error(f"    通道{subchannel.id} {phase}相动量方程求解失败: {e}")
                    v_new = v_old  # 保持原值
                
                # 检查求解结果
                if not np.all(np.isfinite(v_new)):
                    logging.error(f"    通道{subchannel.id} {phase}相求解结果包含无效值")
                    v_new = np.where(np.isfinite(v_new), v_new, v_old)
                
                # 【修复 2025-10-30】使用更温和的速度限制
                # COBRA使用松弛因子而不是硬限制
                # 策略：允许较大变化，但通过松弛因子控制
                
                # 计算相对变化
                dv = v_new - v_old
                v_ref = np.maximum(np.abs(v_old), 0.1)  # 参考速度，最小0.1防除零
                relative_change = np.abs(dv) / v_ref
                
                # 放宽限制：单次迭代最大200%相对变化（原50%太严格）
                max_relative_change = 2.0
                exceed_mask = relative_change > max_relative_change
                if np.any(exceed_mask):
                    dv[exceed_mask] = np.sign(dv[exceed_mask]) * max_relative_change * v_ref[exceed_mask]
                    v_new = v_old + dv
                    if debug_channel:
                        logging.info(f"    {phase}相速度相对变化受限：{np.sum(exceed_mask)}个节点超过200%")
                
                # 【修复 2025-10-30】应用松弛（COBRA风格）
                # 对于稳态问题，COBRA使用阻尼因子DAMPFF
                v_relaxed = (1 - omega) * v_old + omega * v_new
                
                # 将松弛后的速度存入中间场
                self.intermediate_velocity[subchannel.id][phase]['axial'] = v_relaxed
                
                # 【关键修复】同时更新主速度场，否则下次迭代仍用旧值
                subchannel.velocity[phase]['axial'] = v_relaxed.copy()
                
                # 【修复 2025-11-05】使用实际更新后的速度计算残差（修复假收敛问题）
                # 问题：之前用v_new-v_old，但实际更新是v_relaxed-v_old，导致残差虚高
                # COBRA: 直接用实际的新旧场值计算残差
                residual_abs = np.max(np.abs(v_relaxed - v_old))  # 使用松弛后的值
                v_scale = np.max(np.abs(v_old)) + 0.1
                residual = residual_abs / v_scale
                if debug_channel:
                    v_new_residual = np.max(np.abs(v_new - v_old))
                    logging.info(f"    {phase}相动量残差（相对）: {residual:.3e}, 绝对: {residual_abs:.3e}")
                    logging.info(f"    残差诊断 - 理论残差: {v_new_residual:.3e}, 实际更新: {residual_abs:.3e}, 松弛因子: {omega:.3f}")
                
                if phase == 'vapor': max_residual_axial_g = max(max_residual_axial_g, residual_abs)
                else: max_residual_axial_l = max(max_residual_axial_l, residual_abs)

        # 【修复 2025-11-05】修改入口边界条件处理策略（修复假收敛问题）
        # 问题：每次迭代都强制覆盖入口速度，导致入口附近速度永远不变
        # 修复：只在第一次迭代设置，后续让速度自然发展
        # COBRA方法：边界条件在矩阵中处理，不在求解后覆盖
        if hasattr(self, 'boundary_conditions') and 'inlet_velocity' in self.boundary_conditions:
            inlet_bc = self.boundary_conditions['inlet_velocity']
            # NOTE: 只在第一次迭代时强制设置入口边界条件
            if not hasattr(self, '_inlet_bc_iteration_applied'):
                self._inlet_bc_iteration_applied = True
                for subchannel in subchannels:
                    for phase in ['vapor', 'liquid']:
                        if phase in inlet_bc and 'axial' in inlet_bc[phase]:
                            inlet_value = inlet_bc[phase]['axial']
                            # 强制入口节点（j=0）速度为边界条件值
                            subchannel.velocity[phase]['axial'][0] = inlet_value
                            self.intermediate_velocity[subchannel.id][phase]['axial'][0] = inlet_value
                logging.info("  入口速度边界条件已设置（仅第一次迭代）")
                        
        # 【修复 2025-11-04】横向动量方程隐式求解
        logging.info("  开始求解横向动量方程(隐式)")
        lateral_count = 0
        
        # 横向动量方程系数矩阵缓存
        if not hasattr(self, 'momentum_ap_lateral'):
            self.momentum_ap_lateral = {}
            
        for subchannel in subchannels:
            for conn in subchannel.connections:
                if conn['type'] == 'lateral':
                    lateral_count += 1
                    neighbor = next((ch for ch in subchannels if ch.id == conn['neighbor_id']), None)
                    if neighbor is None:
                        continue

                    debug_conn = (subchannel.id == 0 and conn['neighbor_id'] == 1)  # 调试第一个横向连接
                    
                    # 为每个相创建三对角系统
                    for phase in ['vapor', 'liquid']:
                        # 初始化矩阵系数
                        diag = np.zeros(n_axial)
                        rhs = np.zeros(n_axial)
                        lower = np.zeros(n_axial-1)
                        upper = np.zeros(n_axial-1)
                        
                        # 保存旧速度作为参考
                        v_old = subchannel.velocity[phase]['lateral'].copy()
                        
                        if debug_conn:
                            logging.info(f"    调试通道{subchannel.id}-{conn['neighbor_id']} {phase}相横向动量方程:")
                        
                        for i in range(n_axial):
                            # 获取当前节点的物性和空泡份额
                            props = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.VAPOR if phase == 'vapor' else PhaseType.LIQUID)
                            rho = props['density']
                            mu = props['viscosity']
                            alpha = subchannel.alpha[i]
                            v_lat = v_old[i]
                            v_lat_other = subchannel.velocity['liquid' if phase == 'vapor' else 'vapor']['lateral'][i]
                            
                            props_other = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.LIQUID if phase == 'vapor' else PhaseType.VAPOR)
                            rho_other, mu_other = props_other['density'], props_other['viscosity']
                            
                            # 使用小的伪时间步长提供阻尼
                            # 【修复 2025-11-05】进一步降低气相伪时间步长以提高稳定性
                            pseudo_dt = 1e-4 if phase == 'vapor' else 1e-3  # 气相1e-4（从5e-4进一步降低），液相1e-3
                            
                            # --- 时间项 ---
                            time_coeff = 1.0 / pseudo_dt
                            diag[i] += time_coeff
                            rhs[i] += time_coeff * v_old[i]
                            
                            # --- 压力项 ---
                            # 横向压力梯度
                            dp = (neighbor.pressure[i] - subchannel.pressure[i]) / conn['gap_width']
                            if rho > 1e-6:
                                rhs[i] += pseudo_dt * (dp / rho)
                                
                                if debug_conn and i < 3:
                                    logging.info(f"      节点{i}: 压力梯度项={dp:.2e}, 贡献={pseudo_dt*dp/rho:.2e}")
                            
                            # --- 对流项 ---
                            # 简化的对流项处理，只考虑轴向接管
                            if i > 0:
                                v_axial_up = subchannel.velocity[phase]['axial'][i-1]
                                if v_axial_up > 0:  # 上游格式
                                    convection_coeff = v_axial_up / self.grid.dz
                                    lower[i-1] -= convection_coeff
                                    diag[i] += convection_coeff
                            
                            # --- 壁面摩擦项（隐式处理） ---
                            if abs(v_lat) > 1e-6 and rho > 1e-6:
                                # 计算雷诺数
                                Re_lat = rho * abs(v_lat) * conn['gap_width'] / (mu + 1e-10)
                                # 摩擦系数
                                if Re_lat > 2300:  # 湍流
                                    f_lat = 0.316 / (Re_lat**0.25)
                                else:  # 层流
                                    f_lat = 64.0 / (Re_lat + 1e-6)
                                
                                # 【修复 2025-11-05】从配置文件读取横向流动阻力系数
                                lateral_resistance = getattr(self.physics, 'lateral_flow_resistance', 0.7)  # 默认0.7
                                
                                # 隐式摩擦项系数（放在DIAG）
                                friction_coeff = pseudo_dt * f_lat * lateral_resistance * abs(v_lat) / conn['gap_width']
                                diag[i] += friction_coeff
                                
                                if debug_conn and i < 3:
                                    logging.info(f"      节点{i}: 摩擦系数={friction_coeff:.2e}, Re_lat={Re_lat:.0f}, f_lat={f_lat:.4f}")
                            
                            # --- 界面摩擦项（隐式处理） ---
                            if alpha > 1e-6 and (1-alpha) > 1e-6:
                                # 使用物理模型中的界面摩擦系数
                                C_i = self.physics.C_i
                                # 界面面积估算
                                interfacial_area = conn['gap_width'] * self.grid.dz
                                interfacial_coeff = C_i * interfacial_area
                                
                                # 乘以dt统一量纲
                                interfacial_coeff_dt = pseudo_dt * interfacial_coeff
                                
                                if phase == 'vapor':
                                    # 气相：F_i = C_i * (v_g - v_l)
                                    # 隐式：diag += dt*C_i, rhs -= dt*C_i * v_l
                                    diag[i] += interfacial_coeff_dt / (alpha * rho + 1e-10)
                                    rhs[i] += interfacial_coeff_dt * v_lat_other / (alpha * rho + 1e-10)
                                else:
                                    # 液相：F_i = -C_i * (v_g - v_l)
                                    diag[i] += interfacial_coeff_dt / ((1-alpha) * rho + 1e-10)
                                    rhs[i] -= interfacial_coeff_dt * v_lat_other / ((1-alpha) * rho + 1e-10)
                                
                                if debug_conn and i < 3:
                                    logging.info(f"      节点{i}: 界面摩擦系数={interfacial_coeff:.2e}, C_i={C_i:.2f}")

                        # 检查系数矩阵的合理性
                        if debug_conn:
                            logging.info(f"    系数矩阵诊断: diag范围=[{np.min(diag):.2e}, {np.max(diag):.2e}]")
                            logging.info(f"    右端项诊断: rhs范围=[{np.min(rhs):.2e}, {np.max(rhs):.2e}]")
                            # 检查对角占优性
                            diag_dominance = np.min(np.abs(diag)) / (np.max(np.abs(upper)) + np.max(np.abs(lower)) + 1e-10)
                            logging.info(f"    对角占优度: {diag_dominance:.2e} (>1为对角占优)")
                        
                        # 检查矩阵条件数
                        if np.any(np.abs(diag) < 1e-12):
                            logging.warning(f"    通道{subchannel.id}-{conn['neighbor_id']} {phase}相: 系数矩阵对角元过小")
                            diag = np.maximum(np.abs(diag), 1e-10) * np.sign(diag)
                        
                        # 构建三对角矩阵格式
                        ab = np.zeros((3, n_axial))
                        ab[0, 1:] = upper  # 上对角
                        ab[1, :] = diag    # 主对角
                        ab[2, :-1] = lower # 下对角
                        
                        # 缓存系数矩阵
                        conn_key = f"{subchannel.id}_{conn['neighbor_id']}"
                        if conn_key not in self.momentum_ap_lateral:
                            self.momentum_ap_lateral[conn_key] = {'vapor': {}, 'liquid': {}}
                        self.momentum_ap_lateral[conn_key][phase] = diag.copy()
                        
                        # 【修改2025-11-05】应用分离的气液相松弛因子
                        omega = self.relaxation_factors['velocity_g'] if phase == 'vapor' else self.relaxation_factors['velocity_l']
                        
                        try:
                            # 使用三对角求解器
                            v_new = solve_banded((1, 1), ab, rhs)
                            
                            # 检查求解结果
                            if not np.all(np.isfinite(v_new)):
                                logging.error(f"    通道{subchannel.id}-{conn['neighbor_id']} {phase}相求解结果包含无效值")
                                v_new = np.where(np.isfinite(v_new), v_new, v_old)
                            
                            # 应用松弛（COBRA风格）
                            v_relaxed = (1 - omega) * v_old + omega * v_new
                            
                            # 将松弛后的速度存入中间场
                            self.intermediate_velocity[subchannel.id][phase]['lateral'] = v_relaxed.copy()
                            
                            # 同时更新主速度场
                            subchannel.velocity[phase]['lateral'] = v_relaxed.copy()
                            
                            # 【修复 2025-11-05】使用实际更新后的速度计算残差（修复假收敛问题）
                            residual_abs = np.max(np.abs(v_relaxed - v_old))  # 使用松弛后的值
                            v_scale = np.max(np.abs(v_old)) + 0.1
                            residual = residual_abs / v_scale
                            
                            if debug_conn:
                                v_new_residual = np.max(np.abs(v_new - v_old))
                                logging.info(f"    {phase}相横向动量残差（相对）: {residual:.3e}, 绝对: {residual_abs:.3e}")
                                logging.info(f"    残差诊断 - 理论残差: {v_new_residual:.3e}, 实际更新: {residual_abs:.3e}, 松弛因子: {omega:.3f}")
                            
                            if phase == 'vapor': 
                                max_residual_lateral_g = max(max_residual_lateral_g, residual_abs)
                            else: 
                                max_residual_lateral_l = max(max_residual_lateral_l, residual_abs)
                                
                        except Exception as e:
                            logging.error(f"    通道{subchannel.id}-{conn['neighbor_id']} {phase}相横向动量方程求解失败: {e}")
                            # 保持原值
        
        logging.info(f"  处理了{lateral_count}个横向连接")
        return max_residual_axial_g, max_residual_axial_l, max_residual_lateral_g, max_residual_lateral_l

    

  

        # 2. 求解压力修正值 p'
        try:
            p_corr = spsolve(A.tocsr(), b)
            p_corr = p_corr.reshape((n_channels, n_axial))
        except Exception as e:
            logging.error(f"  压力修正方程求解失败: {e}")
            p_corr = np.zeros((n_channels, n_axial))

        # 3. 修正压力和速度
        max_residual = 0.0
        for i in range(n_channels):
            channel = subchannels[i]
            
            # 修正压力
            omega_p = self.relaxation_factors['pressure']
            pressure_old = channel.pressure.copy()
            channel.pressure += omega_p * p_corr[i, :]
            max_residual = max(max_residual, np.max(np.abs(channel.pressure - pressure_old)))

            # 修正轴向速度 (基于中间速度场)
            for j in range(n_axial - 1):
                # 重建 a_P 以计算正确的 d_v
                props_g = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                props_l = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                rho_g, rho_l = props_g['density'], props_l['density']
                rho_mix_j = rho_l * (1 - channel.alpha[j]) + rho_g * channel.alpha[j]
                v_mix_j = self.intermediate_velocity[i]['liquid']['axial'][j]*(1-channel.alpha[j]) + self.intermediate_velocity[i]['vapor']['axial'][j]*channel.alpha[j]
                pseudo_dt = self.pseudo_cfl * dz / max(abs(v_mix_j), 1.0)
                a_p_v = rho_mix_j * channel.area * dz / pseudo_dt
                d_v = channel.area / (a_p_v + 1e-9) # 注意：速度修正不应再乘以松弛因子

                dp_dz = (p_corr[i, j+1] - p_corr[i, j]) / dz
                vel_correction = d_v * dp_dz

                # 两相速度都需要修正
                channel.velocity['vapor']['axial'][j] = self.intermediate_velocity[i]['vapor']['axial'][j] - vel_correction
                channel.velocity['liquid']['axial'][j] = self.intermediate_velocity[i]['liquid']['axial'][j] - vel_correction

            # 修正横向速度 (基于中间速度场)
            for j in range(n_axial):
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        neighbor_id = conn['neighbor_id']
                        # 重建 a_P_lat
                        props_g = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                        props_l = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                        rho_g, rho_l = props_g['density'], props_l['density']
                        rho_mix_j = rho_l * (1 - channel.alpha[j]) + rho_g * channel.alpha[j]
                        v_mix_j = self.intermediate_velocity[i]['liquid']['axial'][j]*(1-channel.alpha[j]) + self.intermediate_velocity[i]['vapor']['axial'][j]*channel.alpha[j]
                        pseudo_dt = self.pseudo_cfl * dz / max(abs(v_mix_j), 1.0)
                        a_p_lat = rho_mix_j * (conn['gap_width'] * dz) * self.grid.dx / pseudo_dt
                        d_v_lat = (conn['gap_width'] * dz) / (a_p_lat + 1e-9)

                        dp_dx = (p_corr[neighbor_id, j] - p_corr[i, j]) / self.grid.dx
                        lat_vel_correction = d_v_lat * dp_dx

                        channel.velocity['vapor']['lateral'][j] = self.intermediate_velocity[i]['vapor']['lateral'][j] - lat_vel_correction
    def _solve_pressure_void_coupled(self, subchannels):
        """
        【PKM简化算法】耦合求解压力和空泡份额
        
        核心改进：参考《Calculation of two-fluid subchannels model of pressurized water 
        reactor: Picard Krylov method》论文，将原PKM-Lite的2×2块矩阵直接求解
        替换为：
        1. 高斯消元降维：2N×2N → N×N（只保留压力方程）
        2. Krylov子空间求解：GMRES + PKM-P物理预条件
        3. 回代求解空泡份额
        
        优势：
        - 维度降低50%，大幅降低系统条件数
        - PKM-P预条件保留轴向强耦合、忽略横向弱耦合，加速收敛
        - GMRES对非对称矩阵更稳定
        
        边界条件处理、速度修正、质量残差计算逻辑保持不变。
        """
        logging.info("  [PKM简化] 开始耦合求解压力和空泡份额（降维+Krylov）")
        n_channels = self.grid.n_channels
        n_axial = self.grid.n_axial
        total_p_unknowns = n_channels * n_axial
        total_unknowns = 2 * total_p_unknowns

        A_pp = lil_matrix((total_p_unknowns, total_p_unknowns))
        A_pa = lil_matrix((total_p_unknowns, total_p_unknowns))
        A_ap = lil_matrix((total_p_unknowns, total_p_unknowns))
        A_aa = lil_matrix((total_p_unknowns, total_p_unknowns))
        b = np.zeros(total_unknowns)

        # 内部节点组装
        for i in range(n_channels):
            for j in range(1, n_axial - 1): # 内部节点
                idx_p = i * n_axial + j  # 压力未知数的索引
                channel = subchannels[i]
                dz = self.grid.dz
                dx = self.grid.dx

                a_p_g_j_minus_1 = self.momentum_ap[i]['vapor']['axial'][j-1]
                a_p_l_j_minus_1 = self.momentum_ap[i]['liquid']['axial'][j-1]
                d_g_axial_j_minus_1 = channel.area / (a_p_g_j_minus_1 + 1e-9)
                d_l_axial_j_minus_1 = channel.area / (a_p_l_j_minus_1 + 1e-9)

                a_p_g_j = self.momentum_ap[i]['vapor']['axial'][j]
                a_p_l_j = self.momentum_ap[i]['liquid']['axial'][j]
                d_g_axial_j = channel.area / (a_p_g_j + 1e-9)
                d_l_axial_j = channel.area / (a_p_l_j + 1e-9)

                d_g_lat, d_l_lat = {}, {}
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        neighbor_id = conn['neighbor_id']
                        props_g = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                        props_l = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                        rho_g, rho_l = props_g['density'], props_l['density']
                        v_mix_j = self.intermediate_velocity[i]['liquid']['axial'][j]*(1-channel.alpha[j]) + self.intermediate_velocity[i]['vapor']['axial'][j]*channel.alpha[j]
                        pseudo_dt = self.pseudo_cfl * dz / max(abs(v_mix_j), 1.0)
                        a_p_g_lat = (1 - channel.alpha[j]) * rho_g * (conn['gap_width'] * dz) * dx / pseudo_dt
                        a_p_l_lat = channel.alpha[j] * rho_l * (conn['gap_width'] * dz) * dx / pseudo_dt
                        d_g_lat[neighbor_id] = (conn['gap_width'] * dz) / (a_p_g_lat + 1e-9)
                        d_l_lat[neighbor_id] = (conn['gap_width'] * dz) / (a_p_l_lat + 1e-9)

                v_g_up = self.intermediate_velocity[i]['vapor']['axial'][j-1]
                v_l_up = self.intermediate_velocity[i]['liquid']['axial'][j-1]
                v_g_down = self.intermediate_velocity[i]['vapor']['axial'][j]
                v_l_down = self.intermediate_velocity[i]['liquid']['axial'][j]

                props_g_up = self.physics.calculate_properties(channel.pressure[j-1], PhaseType.VAPOR)
                props_l_up = self.physics.calculate_properties(channel.pressure[j-1], PhaseType.LIQUID)
                props_g_down = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                props_l_down = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)

                m_dot_g_in_axial = props_g_up['density'] * channel.alpha[j-1] * v_g_up * channel.area
                m_dot_l_in_axial = props_l_up['density'] * (1 - channel.alpha[j-1]) * v_l_up * channel.area
                m_dot_g_out_axial = props_g_down['density'] * channel.alpha[j] * v_g_down * channel.area
                m_dot_l_out_axial = props_l_down['density'] * (1 - channel.alpha[j]) * v_l_down * channel.area

                m_dot_g_in_lat, m_dot_l_in_lat = 0.0, 0.0
                m_dot_g_out_lat, m_dot_l_out_lat = 0.0, 0.0
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        neighbor_id = conn['neighbor_id']
                        v_g_lat = self.intermediate_velocity[i]['vapor']['lateral'][j]
                        v_l_lat = self.intermediate_velocity[i]['liquid']['lateral'][j]
                        props_g_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                        props_l_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                        gap_area = conn['gap_width'] * dz
                        m_dot_g_lat_val = props_g_lat['density'] * channel.alpha[j] * v_g_lat * gap_area
                        m_dot_l_lat_val = props_l_lat['density'] * (1 - channel.alpha[j]) * v_l_lat * gap_area
                        if v_g_lat > 0: m_dot_g_out_lat += m_dot_g_lat_val
                        else: m_dot_g_in_lat -= m_dot_g_lat_val
                        if v_l_lat > 0: m_dot_l_out_lat += m_dot_l_lat_val
                        else: m_dot_l_in_lat -= m_dot_l_lat_val

                R_g_star = (m_dot_g_in_axial - m_dot_g_out_axial) + (m_dot_g_in_lat - m_dot_g_out_lat)
                R_l_star = (m_dot_l_in_axial - m_dot_l_out_axial) + (m_dot_l_in_lat - m_dot_l_out_lat)

                b[idx_p] = R_l_star
                b[idx_p + total_p_unknowns] = R_g_star

                # ===== 【关键修复 2025-10-29】A_aa = ∂R_g/∂α[j] =====
                # 物理意义：气相质量守恒对空泡份额的雅可比
                # 连续性方程：(ρ·α·A)/Δt · Δα + ∂(ρ·α·v·A)/∂z = 0
                # 
                # 【之前的错误】：只有出流项，导致ap_aa全为负数，对角元接近零
                # 【正确组装】：时间项 + 入流项 - 出流项
                
                # 1. 时间项（主导项，保证正定）
                ap_aa = props_g_down['density'] * channel.area / self.grid.dt
                
                # 2. 轴向对流项（迎风格式）
                if j > 0:  # 有上游节点
                    # 入流项（加法）
                    ap_aa += props_g_up['density'] * max(v_g_up, 0.0) * channel.area
                if j < n_axial - 1:  # 有下游节点
                    # 出流项（减法）
                    ap_aa -= props_g_down['density'] * max(v_g_down, 0.0) * channel.area
                
                # 3. 横向流动项
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        v_g_lat = self.intermediate_velocity[i]['vapor']['lateral'][j]
                        gap_area = conn['gap_width'] * dz
                        props_g_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                        # 根据流向决定符号
                        if v_g_lat > 0:  # 流出
                            ap_aa -= props_g_lat['density'] * v_g_lat * gap_area
                        else:  # 流入
                            ap_aa += props_g_lat['density'] * abs(v_g_lat) * gap_area
                
                A_aa[idx_p, idx_p] = ap_aa
                
                # ===== 【诊断2】输出前10个节点的A_aa组装信息 =====
                if idx_p < 10:  # 只输出前10个节点（避免日志过多）
                    logging.info(f"    [诊断] 节点索引={idx_p}(通道{i},轴向{j}): alpha={channel.alpha[j]:.4e}, "
                                f"ap_aa={ap_aa:.4e}, v_g_down={v_g_down:.3f}, v_g_up={v_g_up:.3f}, "
                                f"rho_g={props_g_down['density']:.2f}, P={channel.pressure[j]:.2e}")

                # ===== 【最终修复】A_pa = ∂R_l/∂α[j] =====
                # 液相：ṁ_l = ρ_l·(1-α)·|v|·A
                # ∂R_l/∂α = -∂ṁ_l/∂α = -∂[ρ_l·(1-α)·|v|·A]/∂α = +ρ_l·|v|·A
                
                ap_pa = 0.0
                
                # 轴向通量：使用绝对值
                ap_pa += props_l_down['density'] * abs(v_l_down) * channel.area
                
                # 横向通量：使用绝对值
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        v_l_lat = self.intermediate_velocity[i]['liquid']['lateral'][j]
                        gap_area = conn['gap_width'] * dz
                        props_l_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                        ap_pa += props_l_lat['density'] * abs(v_l_lat) * gap_area
                
                A_pa[idx_p, idx_p] = ap_pa

                # ===== 【关键修复 2025-10-29】采用COBRA的完整组装方法 =====
                # 参考COBRA/BORONC子程序第596行：
                # DIAG = A*DXDT*RHO - COEFA + COEFB + ...
                #      = ρ*A/Δt (时间项，主导) + 入流项 - 出流项
                
                # 物理意义：
                # - 时间项（ρ*A/Δt）提供数值稳定性，保证对角元>0
                # - 入流项（+ρ*v_in*A）：上游流入的质量通量
                # - 出流项（-ρ*v_out*A）：下游流出的质量通量
                # - 速度-压力耦合通过d系数实现
                
                ap_pp, ap_ap = 0.0, 0.0
                
                # 1. 时间项（主导项，参考COBRA）
                # NOTE: 这不是"稳定化"项，而是物理上真实的时间导数项
                ap_pp += props_l_down['density'] * (1 - channel.alpha[j]) * channel.area / self.grid.dt
                ap_ap += props_g_down['density'] * channel.alpha[j] * channel.area / self.grid.dt
                
                # 2. 轴向对流项（迎风格式，参考COBRA第575-588行）
                # 液相
                if v_l_up > 0:  # 正向流：上游流入
                    ap_pp += props_l_up['density'] * (1 - channel.alpha[j-1]) * d_l_axial_j_minus_1 * (1.0 / dz)
                # else: 回流时上游项为零
                
                if v_l_down > 0:  # 正向流：当前节点流出
                    ap_pp -= props_l_down['density'] * (1 - channel.alpha[j]) * d_l_axial_j * (1.0 / dz)
                # else: 回流时下游项在非对角元
                
                # 气相
                if v_g_up > 0:
                    ap_ap += props_g_up['density'] * channel.alpha[j-1] * d_g_axial_j_minus_1 * (1.0 / dz)
                
                if v_g_down > 0:
                    ap_ap -= props_g_down['density'] * channel.alpha[j] * d_g_axial_j * (1.0 / dz)
                
                # 3. 横向流动项（简化处理，参考COBRA DIFFER子程序）
                # NOTE: 常数密度模型，密度效应项d_rho_dp=0，只保留速度耦合
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        neighbor_id = conn['neighbor_id']
                        v_g_lat = self.intermediate_velocity[i]['vapor']['lateral'][j]
                        v_l_lat = self.intermediate_velocity[i]['liquid']['lateral'][j]
                        
                        # 横向速度-压力耦合（迎风格式）
                        if neighbor_id in d_l_lat and v_l_lat > 0:  # 流出
                            ap_pp -= props_l_down['density'] * (1 - channel.alpha[j]) * d_l_lat[neighbor_id] / dx
                        if neighbor_id in d_g_lat and v_g_lat > 0:  # 流出
                            ap_ap -= props_g_down['density'] * channel.alpha[j] * d_g_lat[neighbor_id] / dx

                # 非对角元：上游节点（j-1）耦合（迎风格式）
                # 参考COBRA：COEFC = FJM1*RHOQ(I,JM1) when FJM1 > 0
                if v_l_up > 0:  # 正向流：j-1节点影响j节点
                    A_pp[idx_p, idx_p - 1] = props_l_up['density'] * (1 - channel.alpha[j-1]) * d_l_axial_j_minus_1 / dz
                # else: 回流时此项为零
                
                if v_g_up > 0:
                    A_ap[idx_p, idx_p - 1] = props_g_up['density'] * channel.alpha[j-1] * d_g_axial_j_minus_1 / dz

                # 非对角元：横向耦合（迎风格式）
                # 参考COBRA DIFFER子程序对横向流动的处理
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        neighbor_id = conn['neighbor_id']
                        if self.grid.is_wall_channel(neighbor_id):
                            continue  # 壁面通道无横向耦合
                        
                        v_g_lat = self.intermediate_velocity[i]['vapor']['lateral'][j]
                        v_l_lat = self.intermediate_velocity[i]['liquid']['lateral'][j]
                        
                        # 迎风格式：只在流入时添加非对角元
                        if neighbor_id in d_l_lat and v_l_lat < 0:  # 从邻居流入
                            props_l_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                            A_pp[idx_p, neighbor_id * n_axial + j] = props_l_lat['density'] * (1 - channel.alpha[j]) * d_l_lat[neighbor_id] / dx
                        
                        if neighbor_id in d_g_lat and v_g_lat < 0:  # 从邻居流入
                            props_g_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                            A_ap[idx_p, neighbor_id * n_axial + j] = props_g_lat['density'] * channel.alpha[j] * d_g_lat[neighbor_id] / dx

                A_pp[idx_p, idx_p] = ap_pp
                A_ap[idx_p, idx_p] = ap_ap

        # --- 边界条件处理 ---
        bc_p = self.boundary_conditions.get('pressure_correction', {})
        bc_a = self.boundary_conditions.get('void_fraction_correction', {})

        for i in range(n_channels):
            # 入口 (j=0)
            idx_inlet = i * n_axial + 0
            p_inlet_bc = bc_p.get('inlet', {})
            a_inlet_bc = bc_a.get('inlet', {})

            if p_inlet_bc.get('type') == 'fixedValue':
                A_pp[idx_inlet, idx_inlet] = 1.0
                b[idx_inlet] = p_inlet_bc.get('value', 0.0)
            
            if a_inlet_bc.get('type') == 'fixedValue':
                A_aa[idx_inlet, idx_inlet] = 1.0
                b[idx_inlet + total_p_unknowns] = a_inlet_bc.get('value', 0.0)

            # 出口 (j=n_axial-1)
            idx_outlet = i * n_axial + (n_axial - 1)
            p_outlet_bc = bc_p.get('outlet', {})
            a_outlet_bc = bc_a.get('outlet', {})

            if p_outlet_bc.get('type') == 'fixedValue':
                A_pp[idx_outlet, idx_outlet] = 1.0
                b[idx_outlet] = p_outlet_bc.get('value', 0.0)

            if a_outlet_bc.get('type') == 'zeroGradient':
                idx_prev = idx_outlet - 1
                A_aa[idx_outlet, idx_outlet] = 1.0
                A_aa[idx_outlet, idx_prev] = -1.0
                b[idx_outlet + total_p_unknowns] = 0.0

        # ========================================
        # 【方案A - 2025-10-29】直接求解2×2块矩阵系统
        # ========================================
        # 架构变更：
        # 旧方法（已废弃）：高斯消元降维 → GMRES → 回代
        #   - 要求4×4系统（质量+能量）才能安全降维
        #   - 当前仅2×2系统（无能量方程），A_aa极小导致数值失效
        # 
        # 新方法：直接求解完整2N×2N块矩阵
        #   - 避免A_aa⁻¹的数值问题
        #   - BiCGSTAB+块预条件求解
        # ========================================
        
        logging.info("  [方案A] 采用直接求解2×2块矩阵方法")
        
        try:
            # 直接求解2N×2N系统：[A_pp A_pa; A_ap A_aa] * [Δp; Δα] = [R_l; R_g]
            p_corr_flat, alpha_corr_flat, success = self._solve_2x2_block_system_direct(
                A_pp, A_pa, A_ap, A_aa, 
                b[:total_p_unknowns],  # R_l
                b[total_p_unknowns:],  # R_g
                n_channels, n_axial
            )
            
            if not success:
                logging.error("  [方案A] 2×2块矩阵求解失败")
                return 0.0, 0.0, False
            
            # Reshape为二维数组
            p_corr = p_corr_flat.reshape((n_channels, n_axial))
            alpha_corr = alpha_corr_flat.reshape((n_channels, n_axial))
            
            logging.info("  [方案A] 耦合系统求解成功")

        except Exception as e:
            logging.error(f"  [方案A] 求解过程异常: {e}")
            import traceback
            logging.error(f"  [方案A] 错误堆栈: {traceback.format_exc()}")
            return 0.0, 0.0, False
        
        # ========================================
        # 【已废弃 - 2025-10-29】以下是旧的降维方法（已注释）
        # ========================================
        # logging.info("  [PKM简化] 采用降维+Krylov方法求解压力-空泡耦合系统")
        # try:
        #     # 步骤1: 高斯消元降维（2N×2N → N×N）
        #     A_reduced, b_reduced, A_aa_inv_diag = self._gaussian_elimination_2x2(
        #         A_pp, A_pa, A_ap, A_aa, b[:total_p_unknowns], b[total_p_unknowns:]
        #     )
        #     
        #     # 步骤2: GMRES+PKM-P求解压力方程
        #     p_corr_flat, success_pressure = self._solve_pressure_equation_gmres(
        #         A_reduced, b_reduced, n_channels, n_axial
        #     )
        #     
        #     if not success_pressure:
        #         return 0.0, 0.0, False
        #     
        #     # 步骤3: 回代求解空泡份额
        #     alpha_corr_flat = self._backsubstitute_void_fraction(
        #         p_corr_flat, A_ap.tocsr(), A_aa_inv_diag, b[total_p_unknowns:]
        #     )
        #     
        #     # 步骤4: Reshape
        #     p_corr = p_corr_flat.reshape((n_channels, n_axial))
        #     alpha_corr = alpha_corr_flat.reshape((n_channels, n_axial))
        # except Exception as e:
        #     logging.error(f"  [PKM简化] 求解过程异常: {e}")
        #     return 0.0, 0.0, False
        # ========================================

        # 应用修正
        pressure_residual, void_fraction_residual = 0.0, 0.0
        omega_p = self.relaxation_factors['pressure']
        omega_a = self.relaxation_factors['void_fraction']

        for i in range(n_channels):
            channel = subchannels[i]
            pressure_old = channel.pressure.copy()
            channel.pressure += omega_p * p_corr[i, :]
            pressure_residual = max(pressure_residual, np.max(np.abs(channel.pressure - pressure_old)))

            alpha_old = channel.alpha.copy()
            channel.alpha += omega_a * alpha_corr[i, :]
            channel.alpha = np.clip(channel.alpha, self.min_void, self.max_void)
            void_fraction_residual = max(void_fraction_residual, np.max(np.abs(channel.alpha - alpha_old)))

        logging.info(f"  [PKM简化] 压力修正完成，最大残差: {pressure_residual:.3e}")
        logging.info(f"  [PKM简化] 空泡份额修正完成，最大残差: {void_fraction_residual:.3e}")
        
        # 【诊断】验证压力场是否正确更新
        if n_channels > 0 and n_axial > 3:
            ch0 = subchannels[0]
            logging.info(f"  [诊断] 通道0压力更新验证:")
            logging.info(f"    节点0: 压力={ch0.pressure[0]:.2e} Pa, 修正量={omega_p * p_corr[0, 0]:.3e} Pa")
            logging.info(f"    节点1: 压力={ch0.pressure[1]:.2e} Pa, 修正量={omega_p * p_corr[0, 1]:.3e} Pa")
            logging.info(f"    节点2: 压力={ch0.pressure[2]:.2e} Pa, 修正量={omega_p * p_corr[0, 2]:.3e} Pa")
            # 验证对象引用
            logging.info(f"    对象验证: subchannels[0] is channel: {subchannels[0] is ch0}")

        # TODO: 任务 3.2 - 消除重复代码，重构速度修正逻辑
        # 7. 更新速度场
        omega_u_corr = self.relaxation_factors.get('velocity_correction', 0.1)
        max_dv_dt = 10.0

        for i in range(n_channels):
            channel = subchannels[i]
            p_corr_ch = p_corr[i, :]
            
            # 修正轴向速度
            self._correct_velocity(channel, p_corr, p_corr_ch, 'axial', omega_u_corr, max_dv_dt)

            # 修正横向速度
            for conn in channel.connections:
                if conn['type'] == 'lateral':
                    self._correct_velocity(channel, p_corr, p_corr_ch, 'lateral', omega_u_corr, max_dv_dt, conn=conn)


        # FIXME: 任务 1.2 - 计算并存储真实的质量残差
        mass_imbalance_g = np.zeros((n_channels, n_axial))
        mass_imbalance_l = np.zeros((n_channels, n_axial))

        for i in range(n_channels):
            channel = subchannels[i]
            for j in range(1, n_axial - 1): # 遍历内部节点
                # 计算轴向流入流出
                props_g_up = self.physics.calculate_properties(channel.pressure[j-1], PhaseType.VAPOR)
                props_l_up = self.physics.calculate_properties(channel.pressure[j-1], PhaseType.LIQUID)
                m_in_g_axial = props_g_up['density'] * channel.alpha[j-1] * channel.velocity['vapor']['axial'][j-1] * channel.area
                m_in_l_axial = props_l_up['density'] * (1 - channel.alpha[j-1]) * channel.velocity['liquid']['axial'][j-1] * channel.area

                props_g_down = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                props_l_down = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                m_out_g_axial = props_g_down['density'] * channel.alpha[j] * channel.velocity['vapor']['axial'][j] * channel.area
                m_out_l_axial = props_l_down['density'] * (1 - channel.alpha[j]) * channel.velocity['liquid']['axial'][j] * channel.area

                # 计算横向流入流出
                m_in_g_lat, m_out_g_lat = 0.0, 0.0
                m_in_l_lat, m_out_l_lat = 0.0, 0.0
                for conn in channel.connections:
                    if conn['type'] == 'lateral':
                        neighbor_id = conn['neighbor_id']
                        neighbor = subchannels[neighbor_id]
                        gap_area = conn['gap_width'] * self.grid.dz

                        # 从 i 流向 neighbor
                        v_g_lat_ij = channel.velocity['vapor']['lateral'][j]
                        v_l_lat_ij = channel.velocity['liquid']['lateral'][j]
                        props_g_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                        props_l_lat = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                        
                        m_dot_g_ij = props_g_lat['density'] * channel.alpha[j] * v_g_lat_ij * gap_area
                        m_dot_l_ij = props_l_lat['density'] * (1 - channel.alpha[j]) * v_l_lat_ij * gap_area

                        if v_g_lat_ij > 0: m_out_g_lat += m_dot_g_ij
                        else: m_in_g_lat -= m_dot_g_ij

                        if v_l_lat_ij > 0: m_out_l_lat += m_dot_l_ij
                        else: m_in_l_lat -= m_dot_l_ij

                mass_imbalance_g[i, j] = (m_in_g_axial + m_in_g_lat) - (m_out_g_axial + m_out_g_lat)
                mass_imbalance_l[i, j] = (m_in_l_axial + m_in_l_lat) - (m_out_l_axial + m_out_l_lat)

        # 将残差存储在 self 中，供 solve_steady_state 使用
        self.mass_residual_g = np.abs(mass_imbalance_g)
        self.mass_residual_l = np.abs(mass_imbalance_l)

        return pressure_residual, void_fraction_residual, True # 返回成功状态

    # ========================================
    # PKM简化算法辅助函数（新增）
    # ========================================
    
    
    # ========================================
    # 【方案A - 2025-10-29】直接求解2×2块矩阵系统
    # ========================================
    
    def _solve_2x2_block_system_direct(self, A_pp, A_pa, A_ap, A_aa, R_l, R_g, n_channels, n_axial):
        """
        【方案A】直接求解2×2块矩阵系统（不降维）
        
        架构变更原因：
        - PKM降维方法要求4×4系统（质量+能量），才能安全地消去3个变量
        - 当前程序只有2×2系统（仅质量方程，无能量方程）
        - A_aa极小（~1e-08），降维时A_aa⁻¹会放大误差，导致数值失效
        
        新方法：直接求解完整的2N×2N系统
        
        系统形式：
        ┌         ┐ ┌  Δp  ┐   ┌ R_l ┐
        │ A_pp A_pa│ │      │ = │     │  （液相质量守恒）
        │ A_ap A_aa│ │  Δα  │   │ R_g │  （气相质量守恒）
        └         ┘ └      ┘   └     ┘
        
        参数：
            A_pp, A_pa, A_ap, A_aa: 2×2块矩阵的四个子块（lil_matrix格式）
            R_l: 液相质量残差向量 (N×1)
            R_g: 气相质量残差向量 (N×1)
            n_channels, n_axial: 网格参数
            
        返回：
            p_corr: 压力修正值 (N×1)
            alpha_corr: 空泡份额修正值 (N×1)
            success: 是否成功
        """
        logging.info("    [方案A] 开始直接求解2×2块矩阵系统（2N×2N）")
        
        N = n_channels * n_axial
        
        # 转换为CSR格式（高效矩阵运算）
        A_pp_csr = A_pp.tocsr()
        A_pa_csr = A_pa.tocsr()
        A_ap_csr = A_ap.tocsr()
        A_aa_csr = A_aa.tocsr()
        
        # 【诊断】分析各块矩阵的条件数
        logging.info("    [诊断] 2×2块矩阵条件数分析:")
        for mat_name, mat in [('A_pp', A_pp_csr), ('A_pa', A_pa_csr), 
                               ('A_ap', A_ap_csr), ('A_aa', A_aa_csr)]:
            diag = mat.diagonal()
            diag_nonzero = diag[np.abs(diag) > 1e-15]
            if len(diag_nonzero) > 0:
                cond_est = np.max(np.abs(diag_nonzero)) / np.min(np.abs(diag_nonzero))
                logging.info(f"      {mat_name}: 条件数={cond_est:.2e}, "
                           f"对角元范围=[{np.min(np.abs(diag_nonzero)):.2e}, "
                           f"{np.max(np.abs(diag_nonzero)):.2e}]")
        
        # 构建完整的2N×2N块矩阵
        try:
            # 使用bmat构建块矩阵
            A_full = bmat([[A_pp_csr, A_pa_csr],
                          [A_ap_csr, A_aa_csr]], format='csr')
            
            # 构建右端项
            b_full = np.concatenate([R_l, R_g])
            
            logging.info(f"    [方案A] 系统规模: {A_full.shape[0]}×{A_full.shape[1]}")
            logging.info(f"    [方案A] 非零元: {A_full.nnz}, 稀疏度: {100*(1-A_full.nnz/(A_full.shape[0]**2)):.1f}%")
            logging.info(f"    [方案A] 右端项范数: ||b||={np.linalg.norm(b_full):.3e}")
            
            # 【关键修复】行缩放，平衡各行量级
            diag_A = A_full.diagonal()
            row_scale = 1.0 / np.maximum(np.abs(diag_A), 1e-10)  # 防除零
            row_scale_max = np.max(row_scale)
            row_scale_min = np.min(row_scale)
            
            logging.info(f"    [方案A] 行缩放系数范围: [{row_scale_min:.2e}, {row_scale_max:.2e}]")
            
            # 应用行缩放：D⁻¹·A·x = D⁻¹·b
            D_inv = scipy.sparse.diags(row_scale, format='csr')
            A_full_scaled = D_inv @ A_full
            b_full_scaled = D_inv @ b_full
            
            logging.info(f"    [方案A] 缩放后右端项范数: ||b_scaled||={np.linalg.norm(b_full_scaled):.3e}")
            
        except Exception as e:
            logging.error(f"    [方案A] 块矩阵构建失败: {e}")
            return np.zeros(N), np.zeros(N), False
        
        # 构建块预条件器（对缩放后的矩阵）
        M_inv = self._build_block_preconditioner_scaled(A_full_scaled, n_channels, n_axial)
        
        # 使用GMRES求解（对不定矩阵更鲁棒）
        try:
            logging.info("    [方案A] 开始GMRES求解缩放后的2N×2N系统...")
            
            # GMRES参数
            rtol = 1e-4  # 放宽相对容差
            atol = 1e-8  # 放宽绝对容差
            maxiter = 500  # 增加迭代次数
            restart = 50  # GMRES restart参数
            
            solution, info = gmres(
                A_full_scaled, b_full_scaled,
                M=M_inv,
                rtol=rtol,
                atol=atol,
                maxiter=maxiter,
                restart=restart
            )
            
            if info == 0:
                # 计算缩放前的残差
                residual_scaled = b_full_scaled - A_full_scaled @ solution
                rel_residual = np.linalg.norm(residual_scaled) / (np.linalg.norm(b_full_scaled) + 1e-12)
                
                logging.info(f"    [方案A] GMRES收敛成功！缩放系统相对残差: {rel_residual:.3e}")
                
                # 提取压力和空泡份额修正
                p_corr = solution[:N]
                alpha_corr = solution[N:]
                
                # 诊断修正值范围
                logging.info(f"    [方案A] 压力修正: max={np.max(np.abs(p_corr)):.3e}, "
                           f"mean={np.mean(np.abs(p_corr)):.3e}")
                logging.info(f"    [方案A] 空泡修正: max={np.max(np.abs(alpha_corr)):.3e}, "
                           f"mean={np.mean(np.abs(alpha_corr)):.3e}")
                
                # 验证原始系统的残差
                residual_orig = b_full - A_full @ solution
                rel_residual_orig = np.linalg.norm(residual_orig) / (np.linalg.norm(b_full) + 1e-12)
                logging.info(f"    [方案A] 原始系统相对残差: {rel_residual_orig:.3e}")
                
                return p_corr, alpha_corr, True
                
            elif info > 0:
                logging.warning(f"    [方案A] GMRES未收敛，执行了{info}次迭代")
                
                # 检查部分收敛的结果是否可用
                residual_scaled = b_full_scaled - A_full_scaled @ solution
                rel_residual = np.linalg.norm(residual_scaled) / (np.linalg.norm(b_full_scaled) + 1e-12)
                logging.info(f"    [方案A] 部分收敛残差: {rel_residual:.3e}")
                
                # 如果残差可接受，仍使用结果
                if rel_residual < 0.1:  # 10%相对误差
                    logging.info(f"    [方案A] 残差可接受，继续使用结果")
                    p_corr = solution[:N]
                    alpha_corr = solution[N:]
                    return p_corr, alpha_corr, True
                else:
                    p_corr = solution[:N]
                    alpha_corr = solution[N:]
                    return p_corr, alpha_corr, False
                
            else:
                logging.error(f"    [方案A] GMRES求解失败: info={info}")
                logging.info(f"    [方案A] 尝试回退到直接求解器...")
                # 回退：尝试spsolve
                try:
                    from scipy.sparse.linalg import spsolve
                    solution = spsolve(A_full_scaled, b_full_scaled)
                    logging.info(f"    [方案A] 直接求解器成功")
                    p_corr = solution[:N]
                    alpha_corr = solution[N:]
                    return p_corr, alpha_corr, True
                except Exception as e2:
                    logging.error(f"    [方案A] 直接求解器也失败: {e2}")
                    return np.zeros(N), np.zeros(N), False
                
        except Exception as e:
            logging.error(f"    [方案A] GMRES求解异常: {e}")
            import traceback
            logging.error(f"    错误堆栈: {traceback.format_exc()}")
            return np.zeros(N), np.zeros(N), False

    
    def _build_block_preconditioner_scaled(self, A_scaled, n_channels, n_axial):
        """
        为缩放后的矩阵构建预条件器
        
        策略：简单对角预条件（Jacobi）
        M⁻¹ = diag(A_scaled)⁻¹
        
        参数：
            A_scaled: 缩放后的2N×2N系统矩阵（CSR格式）
            n_channels, n_axial: 网格参数
            
        返回：
            M_inv: 预条件矩阵（对角矩阵）
        """
        logging.info("    [方案A] 为缩放后矩阵构建Jacobi预条件器...")
        
        # 提取对角元
        diag = A_scaled.diagonal()
        
        # 统计对角元
        diag_abs = np.abs(diag)
        diag_min = np.min(diag_abs[diag_abs > 0]) if np.any(diag_abs > 0) else 1e-10
        diag_max = np.max(diag_abs)
        neg_count = np.sum(diag < 0)
        zero_count = np.sum(np.abs(diag) < 1e-15)
        
        logging.info(f"      对角元统计: min={diag_min:.2e}, max={diag_max:.2e}")
        logging.info(f"      负对角元: {neg_count}个, 零元: {zero_count}个")
        
        # 构建预条件器：使用绝对值避免负数问题
        diag_safe = np.where(np.abs(diag) > 1e-10, diag, np.sign(diag) * 1e-10)
        diag_inv = 1.0 / diag_safe
        
        # 返回对角预条件
        M_inv = scipy.sparse.diags(diag_inv, format='csr')
        
        logging.info(f"    [方案A] Jacobi预条件器构建完成")
        
        return M_inv
    
    # NOTE: 以下代码是混乱的，需要手动删除或修复
    # TODO: FIXME - 清理从第1525行开始的混乱代码
    
    def _correct_velocity(self, channel, p_corr, p_corr_ch, direction, omega_u_corr, max_dv_dt, conn=None):
        """
        修正速度场（基于压力修正）
        
        参数：
            channel: 当前通道
            p_corr: 所有通道的压力修正 (n_channels × n_axial)
            p_corr_ch: 当前通道的压力修正 (n_axial,)
            direction: 'axial' 或 'lateral'
            omega_u_corr: 速度修正松弛因子
            max_dv_dt: 最大速度变化率
            conn: 连接信息（横向时使用）
        """
        i = channel.id
        n_axial = self.grid.n_axial
        
        for j in range(n_axial - 1 if direction == 'axial' else n_axial):
            if direction == 'axial':
                # 轴向：dp/dz
                dp = (p_corr_ch[j+1] - p_corr_ch[j]) / self.grid.dz
                a_p_g = self.momentum_ap[i]['vapor']['axial'][j]
                a_p_l = self.momentum_ap[i]['liquid']['axial'][j]
                d_g = channel.area / (a_p_g + 1e-9)
                d_l = channel.area / (a_p_l + 1e-9)
            else:  # lateral
                # 横向：dp/dx
                neighbor_id = conn['neighbor_id']
                dp = (p_corr[neighbor_id, j] - p_corr_ch[j]) / self.grid.dx
                
                # 重新计算横向d系数
                props_g = self.physics.calculate_properties(channel.pressure[j], PhaseType.VAPOR)
                props_l = self.physics.calculate_properties(channel.pressure[j], PhaseType.LIQUID)
                rho_g, rho_l = props_g['density'], props_l['density']
                
                v_mix_j = (self.intermediate_velocity[i]['liquid']['axial'][j] * (1 - channel.alpha[j]) +
                          self.intermediate_velocity[i]['vapor']['axial'][j] * channel.alpha[j])
                pseudo_dt = self.pseudo_cfl * self.grid.dz / max(abs(v_mix_j), 1.0)
                
                a_p_g_lat = (1 - channel.alpha[j]) * rho_g * (conn['gap_width'] * self.grid.dz) * self.grid.dx / pseudo_dt
                a_p_l_lat = channel.alpha[j] * rho_l * (conn['gap_width'] * self.grid.dz) * self.grid.dx / pseudo_dt
                
                d_g = (conn['gap_width'] * self.grid.dz) / (a_p_g_lat + 1e-9)
                d_l = (conn['gap_width'] * self.grid.dz) / (a_p_l_lat + 1e-9)
            
            # 计算速度修正
            vel_corr_g = omega_u_corr * d_g * dp
            vel_corr_l = omega_u_corr * d_l * dp
            
            # 应用修正并限制变化率
            for phase, vel_corr in [('vapor', vel_corr_g), ('liquid', vel_corr_l)]:
                dv = -vel_corr  # 注意符号
                dv_dt = dv / self.grid.dt
                
                if abs(dv_dt) > max_dv_dt:
                    dv = np.sign(dv_dt) * max_dv_dt * self.grid.dt
                
                # 基于中间速度场修正
                channel.velocity[phase][direction][j] = self.intermediate_velocity[i][phase][direction][j] + dv
